Effect of grain dissolution on sloping ground

The static and dynamic stability of natural or constructed slopes can be affected by dissolution or dissolution-like phenomena. Their underlying mechanisms, however, remain unclear. New experimental results and discrete element simulations provide particle-level and macroscale information on the consequences of mineral dissolution on slope behavior. At the microscale, load-carrying grain arches develop around dissolving particles, the porosity increases, and contact force chains evolve to form a honeycomb topology. At the macroscale, while vertical settlements are the prevailing deformation pattern, lateral granular movements that create mass wasting are prominent in sloping ground, even under the quasi-static granular loss. Horizontal grain displacement is maximum at the surface and decreases linearly with the distance from the slope surface to become zero at the bottom boundaries, much like vertical granular displacement along the depth. Sediments with smaller friction angles and steeper slopes experience greater displacement, both vertically and horizontally. Slopes become flatter after dissolution, with the reduction in slope angle directly related to the loss in ground elevation, ΔH/Ho. Yet, because of the porous fabric that results from dissolution, vertical shortening is less than the upper bound, estimated from the loss in the solid mass fraction, ΔH/Ho≈SF. Under water-saturated conditions, the post-dissolution fabric may lead to sudden undrained shear and slope slide.

Test procedure and conditions. The model was slowly flooded with a saturated NaCl brine solution permeated from the bottom (duration: 5 h); brine prevented the dissolution of the salt mixed in with the sand. The gradual dissolution was controlled by progressively lowering the salt concentration in the permeated fluid (duration: 3 days). These slow and progressive changes in salt concentration over the long period ensured all www.nature.com/scientificreports/ soluble particles decreased in size at the same time (a "homogeneous dissolution" as in Cha and Santamarina 47 ). The low reaction rate (low Damkohler number) guaranteed uniform dissolution and avoided preferential dissolution modes such as localized dissolution and face dissolution initiated from the inlet boundary. In addition, a very low hydraulic gradient (i < 0.1) was maintained for the saturation and flow so that it created negligible uplift pressure on the sediment particles. A camera recorded the slope deformation by high-resolution timelapse imaging with a five-minute time interval. The experimental setup was isolated on a heavy, vibration-free table during dissolution. While it is impossible to simulate the time scale of dissolution in the field, we ensured that the laboratory dissolution was very gradual, as demonstrated by the long and uniform dissolution. Consequently, the system behaved quasi-statically.
Previous research has used the inertial number (I) to quantify the dynamic effects in particulate discrete element simulations 48 and experimental analysis 49 . The calculated inertial number was ~ 10 −8 for the particles on the surface, which is well within quasi-static criteria (I < 10 −3 ).
Experimental results. Digital image correlations failed to provide adequate displacement fields-from the start to the end-due to a lack of contrast in the sediment. Instead, we manually traced displacements of clearly visible and persistent impurities detectable in high-resolution images. The results in Fig. 2a show the displacement vectors from the beginning to the end of dissolution, which demonstrate that vertical displacements prevail, the magnitude of vertical displacements increases almost linearly with elevation, and horizontal displacements become increasingly prominent towards the slope surface. An areal analysis of the specimen shows that the total volume decreases by 7.5%, with the dissolved fraction of 12% by solid volume (10% by weight), which provides the calculation of the increase in porosity from the initial porosity n 0 = 0.412 to final porosity n f = 0.440.
Notably, we detected a sudden displacement when all the recorded images were played as a movie sequence. Figure 2b presents the two consecutive images that highlight this sudden displacement. Digital image correlation did capture the sudden movement between these consecutive images; see Fig. S1 in the Supplemental Data. The displacement vectors defined a shallow failure surface. The maximum possible duration for this local failure is the time interval between photographs Δt = 5 min (although the fact that the localization is solely contained within the two consecutive images indicates the actual duration could be arbitrarily much shorter). Given a drainage length similar to the slide depth, d = 10 cm, the coefficient of consolidation for drained conditions should be higher than c v > d 2 /Δt = 0.3 cm 2 /s, which corresponds to fine sands. These results suggest that this local failure took place under undrained conditions.
A possible mechanism leading to the slide could involve the following sequence of events, as also illustrated in Fig. 3. Granular dissolution results in increased porosity, as shown in this experiment and other work, and the sediment becomes contractive 47 . At the same time, internal stress states approach the Coulomb failure condition, i.e., K a 50 . Therefore, the active failure condition initiated the structural collapse of the porous granular skeleton. This sudden collapse combined with the highly contractive tendency of the post-dissolution sediments 47 can lead to the temporal generation of excess pore pressure, as described in the above paragraph, and a decrease in the effective stress, which lowered shear strength and resulted in the slope failure/slide.
We used low-friction, scratch-free acrylic plates for the lateral walls. Nevertheless, the boundaries may create friction between the particles and lateral walls and could lessen granular movement at the boundaries to some degree. First, during dissolution, friction may reduce settlement and then once shear localization initiates, it will potentially reduce the amount of movement. Indeed, our shear slide occurred over a short distance. The fact that the localization occurred in spite of the boundary friction may indicate a more substantial slide in a 3D experiment with a large lateral dimension.

Discrete element simulations
The consequences of dissolution on sloping ground were studied further using two-dimensional (2D) discrete element method (DEM) simulations via Particle Flow Code (PFC) 2D by Itasca. Unlike the experiments, the DEM simulations were dry or drained (no fluid module). These conditions still have relevance in dissolution under unsaturated conditions (e.g., by repeated rainwater infiltrations) 51,52 , and solifluction in partially saturated ground, or even saturated slopes of coarse-grained sediments that do not experience localization during dissolution. While not emulating the exact conditions, the simulations and experiments gather complementary information.
Simulation environment and geometry. Table 1 shows the basic simulation environment and material properties. The 2D simulations used the linear contact model with normal stiffness k n = 10 8 N/m and shear stiffness k s = 10 8 N/m. As noted above, conditions were dry/drained, and there was no side friction against the front or back walls. Additionally, there is no friction between the sediment and the vertical lateral walls, which allowed the 1D settlement of the level ground portion subjected to uniform dissolution. The friction coefficient between the sediment and the lower boundary (non-dissolving base material) was 0.5 (Fig. 4f). The stiffnesses of the contact model and the inter-particle friction were selected based on the literature and the software manual [53][54][55][56] . The sediment consisted of disks that were placed at random locations within the model geometry and allowed to grow to their final target size (i.e., uniform size distribution with R min = 0.4 mm and R max = 0.6 mm).
Grain angularity and interlocking, which imply rotational resistance, can be numerically implemented by imposing rolling resistance [57][58][59][60] or by the introduction of non-spherical particles and granular clusters 61-64 . Hindered particle rotation is a computationally efficient approach that accounts for grain angularity 65,66 . Numerical results are physically inconsistent when all particles have hindered rotation; in this study, we hindered the rotation of a preselected percentage of randomly located particles, that is, hinder rotation (HR) of 0%, 40%, and 80%,  47 . Although this approach did mimic the macroscopic trends and global-scale micromechanical parameters, it created more unnatural behaviors than those from computationally expensive techniques that use realistic particle shapes; in particular, we observed excessive, localized dilation around rotationally hindered particles. Simulations used three initial slope angles: β = 20°, 30°, and 40° (Fig. 4). The numerically measured angle of repose, Φ r , increased with the percentage of particles with hindered rotation: Φ r = 22° for HR = 0%, Φ r = 35° for HR = 40%, and Φ r = 48° for HR = 80% 67 . Therefore, only the slope with β o = 20° was stable when HR = 0%; the slopes with β o = 20° and 30° were stable for HR = 40%; and the three slopes with β o = 20°, 30°, and 40° were stable when HR = 80%. The specimen dimensions for the β = 20° cases are a base length of 284 mm and height of 72 mm. Dissolution simulation. Soluble particles were randomly selected and, thus, homogeneously distributed in the preformed granular packing (soluble particles shown in red; Fig. 4). The mass fraction of soluble particles (SF) accounted for SF = 25% of all particles. Note soluble fractions in a given sediment that contains watersoluble minerals such as carbonate and evaporite can vary greatly 68,69 . www.nature.com/scientificreports/ Smooth and gradual size reduction prevented numerical instability and dynamic effects. The dissolution of soluble particles was performed by gradually and simultaneously reducing the radius of all the soluble particles at the same rate 70 . In particular, the gradual size reduction involved numerous steps of minute size reduction (specifically, a radius reduction of 1/50,000 times the initial radius in each step), with each step followed by a full equilibrium stage 47,71 . In addition, the ratio of the mean unbalanced force to the mean contact force was always smaller than 0.001, which ensured stable conditions throughout the dissolution process.
The inertial number I is the ratio between the time for a given displacement when accelerated by the stressdependent skeletal forces σ' d 2 and the time for the same displacement given an imposed strain rate γ 48,72 : Stress ratio K    www.nature.com/scientificreports/ For particles of diameter d = 1 mm, grain density ρ = 2650 kg/m 3 , average effective stress σ' = 0.02 kPa of surface particles (the most critical case), and particle shrinking rate γ = 0.017/s, the computed inertial number I≈2 × 10 −4 is within the quasi-static criterion I < 10 −3 for strain rate independent frictional resistance. For the global movement, applying a slope movement rate γ = 0.004/s and an effective stress 0.6 kPa (at mid depth), the calculated I is 8 × 10 −6 . Macroscale parameters stabilized when the sizes of dissolvable particles decreased to 10% of the initial sizes. We ended our simulations when the particle sizes reduced to 1% of their initial sizes. The time step was set to �t = 0.6 √ m/K , where m is the particle mass and K the particle stiffness; the time step changed as the grain mass m diminished during the simulations. The physical time for the dissolution simulations was about 2 min while the computation time exceeded several weeks. Observations from the simulation results follow.
Contact force and internal fabric. Contact forces form characteristic "chains" in granular materials.
Force chains are evenly distributed before dissolution in our simulations and have a preferentially vertical orientation in agreement with the typical stress ratio K 0 ≈ 0.45 under normal compaction, yet the principal direction of force chains deviate from the vertical near the slope (Fig. 4). Figure 5 shows that higher stress anisotropy exists in the sloping ground (from the surface to the base) in its initial condition. The larger slope angle creates higher values of maximum shear stresses divided by the mean normal stresses.
Force redistribution starts after the dissolvable particles begin to contract; forces initially carried by dissolvable grains transfer to neighboring grains during dissolution (see the video-Supplementary Data). Preponderant interparticle force chains form arches around developing voids in a honeycomb-like fabric (Fig. 6). The honeycomb-shaped network of force chains was more prominent in sediments with more granular interlocking, that is, packings with a higher fraction of rotationally hindered particles, HR (Fig. 6). The coordination number decreased for all cases: from an initial value of 3.62 (all cases), to 2.38 (no hindered rotation), 2.28 (HR = 40%) and 2.09 (HR = 80%) at the end of dissolution. The state of stress became more stable after dissolution, especially for the mid slope. Careful inspection of the images in Fig. 6 reveals large remnant voids located next to major contact force chains after dissolution; this implies that force arches developed around the dissolving particles. A 2D crosscorrelation analysis between the image of large voids ( Fig. 7a; Note: Grains are grown ΔR/R≈50% without displacement to close all small voids) and an image of contact force chains (Fig. 7b) confirms that most large voids were found one particle diameter away from the major force chains (and typically below; see Fig. 7c).
This distinct post-dissolution internal fabric anticipates a different sediment response upon further loading, e.g., more prone to buckling under shear perturbation. In particular, the inherent shear loading in a sloping ground may trigger static liquefaction given the high contractive tendency in soils that experienced dissolution 47 . This may explain the shallow slide reported in Fig. 2b. Displacement. In agreement with the experimental results (Fig. 2a), vertical settlement was the prevailing global deformation pattern (Fig. S2 in the Supplemental Data). Sediments with higher soluble fractions and lower interlocking or macroscopic friction angles experienced larger amounts of settlement. Once again, vertical settlements increased quasilinearly toward the surface because of the random distribution of dissolvable grains. The third column in Table 2 shows vertical strains determined by surface settlements near the vertical walls (i.e., away from the slope). All vertical strains are well below 0.25, in agreement with increased porosity. Higher granular interlocking reduces vertical strain. In addition, vertical strains increase with the initial heights by a small degree due to gains in self-weight with increasing depth ( Table 2).
The sloping ground also creates significant mass wasting, even under the fully quasi-static dissolution. As the randomly distributed soluble particles dissolve and surrounding particles are rearranged, all of the particles are mobilized vertically and horizontally (see the video-Supplementary Data). Their mobilizations are coupled during dissolution in a way that as much as particles displace vertically, horizontal movements co-occur in the sloping ground, which bears large maximum shear stresses relative to mean normal stresses (Fig. 5). Figure 8 presents the horizontal components extracted from the final displacement vectors. Horizontal displacements were largest on the slope surface and decreased away from the slope surface without a sharp transition, that is, no shear localization. The horizontal displacements diminish gradually toward the level ground. Figure 9 quantifies the horizontal granular displacements versus the initial distance perpendicular to the slope surface in the mid slope. In all cases, horizontal grain displacement is maximum at the surface and decreases quasi-linearly with the distance from the slope surface, and becomes zero at the bottom boundaries, similar to the vertical granular displacement along the depth. It resembles simple shear flow although their movements are quasi-static.
Granular interlocking reduces horizontal displacement, and the larger initial slope angle promotes greater mass wasting (Figs. 8 and 9). Clearly horizontal displacements increase with the slope angles and decrease with www.nature.com/scientificreports/ granular interlocking more sensitively than the vertical displacements (Table 2 and Fig. 10). In addition, granular dissolution in sloping grounds increases vertical strain more than dissolution on the level ground (Table 2: compare the third and fourth columns); as grains displace horizontally, the slope geometry may create additional downward movement. www.nature.com/scientificreports/ Histograms of horizontal displacements after dissolution clearly indicate that the bulk of the medium had a mean horizontal displacement near zero (Fig. 11). The green lines show histograms for the ~ 2000 particles (19% of the total number of particles) in the level ground, closest to the rigid wall on the right. In this zone, the horizontal displacements of particles are small, contained within ± 2 mm without a preferential direction; the histogram is nearly symmetric about zero displacement (Fig. 11). When all of the particles are plotted, however, significant portions are displaced horizontally outward (positive direction) (Fig. 11), which represents the grain   www.nature.com/scientificreports/

Discussions
Field implication: mass wasting due to granular mass loss. Certainly, the slope geometry and soluble distribution in this DEM study represent a particular case: a young sloping deposit with a limited slope extent such as energy waste or embankment that contains a soluble fraction overlies an insoluble base ground. Some field situations show that the slope surface and the lower end of the dissolution front may be parallel and the affected depth determined by the intensity and duration of each infiltration event. Rainwater infiltration may create secondary porosity in a young sloping deposit and weaken cementation in sloping sediments. In mass wasting events induced by freezing-thaw, the freezing and thawing depths are determined by the temperature difference. Layer(s) can also dominate a dissolution depth or geometry. Nevertheless, this study shows that horizontal/vertical grain movements are proportional to the affected depth of the dissolution. Also, the effects of major parameters are qualitatively recognized from this study and previous studies. While further, focused studies are required to accurately determine each effect to be useful for predicting the behavior, the following relation lends perspective to the combined effects in the absence of shear localization: where each function correlates positively with the parameter inside the parentheses. θ is the slope angle, SF the fraction of randomly distributed soluble particles, ϕ the macroscopic friction angle, L the dimension of a slope, either the height or length, and D (%) the progression of dissolution that can range from 0 to ~ 80%, after which the movement begins to stabilize.
Final slope angles. The initial slope angle, β 0 , became flatter as dissolution progressed, and the final slope was β final < β 0 in all cases (Fig. 12). In the absence of undrained failure, lateral displacements during dissolution had a minor effect on the slope angle; in fact, a good approximation to the final slope angle was trigonometrically computed from the shortening, ΔH, of the initial land elevation, H o : Apparently, the proximity of the initial slope angle, β 0 , to the numerically measured angle of repose, β repose , has a secondary effect; however, granular interlocking does reduce the vertical shortening (i.e., smaller ΔH/H o ) and preserves steeper slopes after dissolution than less interlocked specimens, as predicted by the equation above. The soluble fraction, SF, is an upper bound estimate for vertical shortening, ΔH/H o < SF-in fact, ΔH/H o = SF when dissolution takes place at a constant porosity (dotted line in Fig. 12); otherwise, porosity increases during dissolution, as observed in this study.
The limited length of slope in this DEM study, as is often the case for embankments and engineered fills, creates a non-uniform slope movement. Horizontal granular displacement is more restricted near the base than on the mid slope because of the non-dissolving base material and its friction with grains of the sloping deposit. The horizontal displacement also diminishes gradually toward the level ground. Non-uniform horizontal movements thus result in slightly curved surfaces (more pronounced in high-angle slopes, e.g., Figs. 6f and 8f).

Conclusions
Mineral dissolution and precipitation are concurrent soil processes and can take place in relatively short time scales in advective regimes and when minerals are exposed far from equilibrium, which is common in young natural and engineered systems. We investigated the effect of grain dissolution on slope stability using a 1 g model experiment complemented by discrete element simulations.
Upon dissolution, load-carrying grain arches develop around dissolving particles. A honeycomb-shaped contact force chain topology characterizes the high-porosity fabric after dissolution. The cross-correlation analysis confirms that large voids remain next to major contact force chains after dissolution. The sloping ground has high initial values of maximum shear stress relative to mean normal stress. After dissolution, the state of stress became more stable. www.nature.com/scientificreports/ Dissolution leads to significant slope movements. While global vertical settlement is the prevailing deformation pattern, grains displace horizontally in the sloping ground, which leads to significant mass wasting even though the dissolution takes places quasi-statically. Horizontal granular movements are maximum at the slope surface and decrease linearly away from the slope surface bounded by the non-dissolving base, similar to vertical granular displacement along the depth. Sediments with a lower global friction angle or steeper slope experience larger displacements, both vertically and horizontally. Distance from the surface [mm] Horizontal displacement [mm] (c) β=40°B oundary Figure 9. The horizontal granular displacements versus the initial distance perpendicular to the slope surface in the mid-slope. The green strip (4 mm width) shows the range of selected grains.  Fig. 9) normalized by the initial heights from the base. Numerical data from www.nature.com/scientificreports/ Slopes become flatter after dissolution; the reduction in slope angles directly related to the loss in ground elevation, ΔH/H o . Yet, for all cases, the slope angle is higher than the estimate from the upper bound for vertical shortening, ΔH/H o ≈SF, because of increased porosities.
No shear localization or catastrophic failure was observed in the numerical simulations during dissolution under drained or dry conditions. However, sudden undrained shear may take place and lead to shallow failures. The sequence of events involves: (1) dissolution and increase in porosity, (2) internal stress approaching failure, (3) grain skeleton collapse, (4) pore pressure buildup, and (5) shear occurring as a form of static liquefaction. www.nature.com/scientificreports/

Data availability
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.